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The local structure of liquid water as a function of temperature is a source of intense research. This structure 
is intimately linked to the dynamics of water molecules, which can be measured using Raman and infrared 
spectroscopies. The assignment of spectral peaks depends on whether they are collective modes or single 
molecule motions. Vibrational modes in liquids are usually considered to be associated to the motions of 
single molecules or small clusters. Using molecular dynamics simulations we find dispersive optical phonon¬ 
like modes in the librational and OH stretching bands. We argue that on subpicosecond time scales these 
modes propagate through water’s hydrogen bond network over distances of up to two nanometers. In the 
long wavelength limit these optical modes exhibit longitudinal-transverse splitting, indicating the presence of 
coherent long range dipole-dipole interactions, as in ice. Our results indicate the dynamics of liquid water 
have more similarities to ice than previously thought. 


The local structure and dynamics of liquid water as a 
function of temperature remains a source of intense re¬ 
search and lively debate>i“— A thus far unrecognised dis¬ 
crepancy exits between the peak assignments reported 
in Raman spectra with those reported in dielectric/IR 
spectra. Although early experimentalists fit the Ra¬ 
man librational band with two peaks^^ it is better fit 
with three (Supplementary Table Previously these 

three peaks were assigned to the three librational motions 
of the water molecule - twisting (ai 435 cm“^), rocking 
(~ 600 cm“^) and wagging (?« 770 cm“^)i2ii2iiii However, 
when comparing these assignments to infrared and di¬ 
electric spectra, one runs into a serious discrepancy. One 
expects to find the two higher frequency modes to be 
present, since only the rocking and wagging librations 
are IR active. The twisting libration, consisting of a ro¬ 
tation of the hydrogen atoms around the C2 axis, is not 
IR active since it does not affect the dipole moment of 
the molecule. Instead, IR spectra show two peaks at 380 
and 665 cm“^fiS and similarly dielectric spectra show 
peaks at 420 and 620 cm“^f.I^ in disagreement with this 
assignment. 

The assignment of longitudinal optical phonon modes 
to Raman spectra can be made by looking at the longi¬ 
tudinal dielectric susceptibility. This method has been 
used previously to assign longitudinal phonon modes to 
the Raman spectra of ice Ih,^®^— ice Ic^ and vitreous 
Ge02 and Si02'^ It has previously been shown that the 
librational peak in the longitudinal dielectric suscepti¬ 
bility of water is dispersive;^ and Bopp & Kornyshev 
noted that the dispersion relation has the appearance 
of an optical phonon mode^ The longitudinal mode in 
the dielectric susceptibility is equivalent to the dispersive 
mode discovered by Ricci et al. (1989) in the spectrum 


of hydrogen density fluctuations 

Comparison of peak positions in longitudinal and 
transverse dielectric susceptibilities often reveals 
longitudinal-transverse (LO-TO) splitting. LO-TO split¬ 
ting indicates the presence of long-range dipole-dipole 
interactions in the system. One way to understand 
LO-TO splitting is through the Lyddane-Sachs-Teller 
(LST) relation)^ 


4 ^=— ( 1 ) 
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Although this relation was originally derived for a cu¬ 
bic ionic crystal it was later shown to have very general 
applicability ) and has been applied to disordered and 
glassy solids.— To apply this equation to water we 
must use a generalized LST relation which takes into ac¬ 
count all of the optically active modes in the system and 
the effects of dampening^2.4 The generalized LST relation 
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Here the index i runs over the Debye peaks in the sys¬ 
tem and the index j runs over the number of damped 
harmonic oscillator peaks. The longitudinal frequencies 
of the damped harmonic oscillators must be considered 
as complex numbers (oJLi = <^Li + * 74 ), where 7 i is the 
dampening factor. 

As shown by Barker, the generalized LST equation 
can be understood purely from a macroscopic point of 
viewj^ so by itself it yields little insight into microscopic 
dynamics. LO-TO splitting can be understood from a 
microscopic standpoint via the equation— 
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FIG. 1. Dielectric susceptibilities of ice and water. 

Computed from index of refraction data using equations U 
and[6l data from 210 to 280 K comes from aerosol droplets^ 
while the data at 300 comes from bulk liquid. 



FIG. 2. Polarization correlation functions. Longitudinal 
(a) and transverse (b) polarization correlation functions (see 
equation [ 3 ) for TIP4P/£, a rigid model. The oscillations at 
small k come from the collective librational mode, which is 
much more pronounced in the longitudinal case. 


II. LO-TO SPLITTING FROM EXPERIMENTAL DATA 


Here v is the volume per unit cell, Qk is the normal 
coordinate of mode fc, and C is a prefactor which de¬ 
pends on the type of lattice and the boundary conditions 
of the region being considered (for an infinite cubic lat¬ 
tice, (7=1). Equation [3] shows that LO-TO splitting is 
intimately related to crystal structure, and it has been 
used to evaluate the quasi-symmetry of room tempera¬ 
ture ionic liquids)^. 

In this work we show how the dielectric susceptibil¬ 
ity can be used to probe water’s local structure and dy¬ 
namics. Our work solves the aforementioned peak as¬ 
signment discrepancy. We find that the lowest frequency 
librational Raman peak {fv 435 cm“^) is a transverse op¬ 
tical phonon-like mode while the highest frequency peak 
(~ 770 cm“^) is a longitudinal optical phonon-like mode. 
This explains why the highest frequency Raman mode 
does not appear in IR or dielectric experiments, since 
such experiments only report the transverse response. 
We show that the transverse counterpart also exhibits 
dispersion. We argue that these dispersive modes are 
due to optical phonons that travel along the H-bond net¬ 
work of water. Our results indicate that not only does 
water exhibit LO-TO splitting, but also that its depen¬ 
dence with temperature is anomalous. We suggest that 
this measurement provides an alternative probe to eval¬ 
uate structural changes in liquid water as a function of 
temperature. 


We wish to study the k dependence of the dielectric 
susceptibility, where k = ‘I'KjX. k— dependence cannot 
be probed directly by experiment, but in the limit of 
infinite wavelength {k —> 0) the longitudinal and trans¬ 
verse dielectric susceptibilities can be obtained from the 
dielectric function via the following relations 


XL{k -k 0,w) = 1 

e{uj) 

(4) 

XT{k -k 0, w) = e{uj) - 1 

(5) 


Note that the transverse susceptibility is what one nor¬ 
mally calls susceptibility. The dielectric function can be 
obtained from the index of refraction n(a;) and extinction 
coefficient k{u}) as: 

e'{u})=n^{uj)-k^{uj) 
e"(a;) = 2n{oj)k{ui) 

These equations allow us to use previously published ex¬ 
perimental dat a^^i^^ to calculate the imaginary part of 
the longitudinal response. We find significant LO-TO 
splitting in the librational and stretching bands (fig. [!])• 

A. Polarization correlation functions 


I. RESULTS 


As in our previous work^ we compared results from a 
rigid (TIP4P/e) model, a flexible model (TIP4P/2005f), 
and a flexible and polarizable model (TTM3F) in all of 
our analyses. 


The normalized longitudinal and transverse polariza¬ 
tion correlation functions are defined as: 

_ (Pi/T(fc,t)-Pl/r(fc,0)) 

- {PL/T{k,0)-Pl/^{kM 

The correlation functions found for TIP4P/e at small 
small k are shown in figure O Since TIP4P/e is a rigid 
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FIG. 3. Imaginary part of longitudinal (a) & trans¬ 
verse (top) susceptibility From a simulation with TTM3F 
at 300 K. In the longitudinal spectra both the librational (« 
750 cm“^) and OH stretching peak (~ 3500 cm“^) peaks ex¬ 
hibit dispersion. 
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TABLE I. Resonance frequencies and lifetimes Frequen¬ 
cies are given in cm“^ and lifetimes in picoseconds. The val¬ 
ues from simulation were computed at the smallest k in the 
system. The experimental values are based on the position of 
the max of the band and therefore only approximate. 


model, only librational motions are present. The addi¬ 
tion of flexibility and polarizability add additional high 
frequency oscillations to the picture (Supplementary Fig. 
1). In the small wavenumber regime (k < 1.75A) there is 
a damped oscillation which corresponds to the collective 
librational phonon-like mode. This damped oscillation is 
superimposed on an underlying exponential relaxation in 
both the transverse and longitudinal cases. In the lon¬ 
gitudinal case the relaxation time t(/c) of the underlying 
exponential relaxation exhibits non-monotonic behaviour 
with k, reaching a maximum at fc ~ 3A (Supplemen¬ 
tary Fig. 2). At wavenumbers greater than k ~ 2.5A 
only intramolecular motions contribute. 


B. Dispersion of the librational peak 

Figure [3] shows the imaginary part of the longitudinal 
and transverse susceptibility for TTM3F. In the longi¬ 
tudinal case the librational peak is clearly seen to shift 
with k. In the transverse case, the lower frequency por¬ 
tion of the band is seen to shift slightly with k. Dis¬ 
persion relations for the longitudinal and transverse li¬ 
brational peaks are shown in figure [S] for three different 
temperatures, using one peak fits. The dispersion rela¬ 
tions appear to be that of optical phonons. In both the 
longitudinal and transverse case the dampening factors 
remain less than the resonance frequencies, indicating an 
underdamped oscillation (Supplementary Fig. 3). The 
longitudinal dispersion relation for TIP4P/2005f agrees 
with that found by Bopp & Kornyshev (who used the 
flexible BJH model).Resat et al. also obtained a simi¬ 
lar dispersion relation (but at a higher frequency), using 
the reference memory function approximation for TIP4P 
instead of molecular dynamics 

Resonance frequencies and lifetimes for the smallest 
k are shown in table III Bl The speed of propagation of 
these modes was computed by finding the slope dio/dk 
in the regime of linear dispersion. For TIP4P/2005f we 
found speeds of ~ 2700 m/s and ~ 1800 m/s for the 
longitudinal and transverse modes. These propagation 
speeds are above the speed of sound in water (1500 m/s) 
but below the speed of sound in ice (4000 m/s). The 
temperature dependence of the propagation speed was 
found to be very small. 

In both the longitudinal and transverse cases, the 
residual of the peak fitting shows features not captured 
by our Debye + one resonance fit of the librational peak. 
In both the longitudinal and transverse cases there is a 
non-dispersive peak at higher frequency, located at ~900 
cm-i for TIP4P/2005f and at a:650 cm"! in TTM3F. 
This peak is negligibly small in the k = 0 longitudinal 
susceptibility but appears as a shoulder as k increases. 
In the transverse case the overlapping peak persists at 
fc = 0, so we found that the k = 0 transverse spectra is 
best fit with two peaks, in agreement with experimen¬ 
tal spectra. As we describe later, the higher frequency 
transverse peak is largely due to the self part of the re¬ 
sponse and is associated with the wagging librations of 
single molecules. 

C. Importance of polarizability 

There are several notable differences between TTM3F 
and the non-polarizable model TIP4P/2005f. First of 
all, the librational band of TIP4P/2005f is at higher fre¬ 
quency, in worse agreement with experiment. This dif¬ 
ference in frequency is likely related to the parameters 
of TIP4P/2005f and not its lack of polarization. More 
importantly, we find that TTM3F exhibits dispersion in 
the OH stretching band (~ 3500 cm“^) in the longitu¬ 
dinal case while TIP4P/2005f does not. The transverse 












FIG. 4. Imaginary part of the longitudinal suscepti¬ 
bility. For TIP4P/2005f at 300 K. No significant dispersion 
is observed in the OH stretching peak. 


susceptibility of TTM3F does not exhibit such dispersion 
but the magnitude of the OH stretching band increases 
at small fc, indicating long range intermolecular correla¬ 
tions. TIP4P/2005f does not exhibit this behavior. Sim¬ 
ilarly, at fc = 0 TTM3F exhibits significant LO-TO split¬ 
ting in the OH stretching band while TIP4P/2005f does 
not (fig. in]). These findings are consistent with Hey- 
den et al.’s results for the A:-resolved IR spectra from ab- 
initio simulation, where they concluded that polarization 
allows for intermolecular correlations at the OH-stretch 
frequency!^ 

These findings can be understood from the dipole 
derivative in equation [3| In the librational band the 
derivative of the dipole moment with respect to normal 
coordinate is purely due to rotation, while in the OH- 
stretching band it is due to changes in the geometry of the 
molecule and electronic polarization of the molecule dur¬ 
ing the OH stretching. In principle there may be coupling 
between the librational and stretching motions, but typi¬ 
cally such rotational-vibrational coupling effects are neg¬ 
ligibly smallj^ The dipole moment surface (fluctuating 
charges) and polarization dipole incorporated in TTM3F 
account for the changes in polarization that occur during 
OH stretching motion. These results confirm the sig¬ 
nificance of polarization in capturing the OH stretching 
response of water 

Figure m shows a comparison of TTM3F, TIP4P/2005f 
and experiment at k = 0. While the location of the peaks 
in TTM3F are in good agreement with the experimen¬ 
tal data at 298 K, the magnitude of the longitudinal re¬ 
sponse is greatly overestimated in TTM3F. The degree of 
LO-TO splitting in the OH stretching peak is also over¬ 
estimated in TTM3F. In general it appears that TTM3F 
overestimates the dipole derivative in equation [3] while 
TIP4P/2005f underestimates it. Figured also shows the 
effect of polarization at low frequencies, in particular the 
appearance of an H-bond stretching response at ~ 250 
cm“^ in TTM3F which is absent in TIP4P/2005fi^ 



FIG. 5. Dispersion relations for the propagating libra¬ 
tional modes. For TIP4P/2005f at three different tempera¬ 
tures (squares = longutudinal, pluses = transverse). A similar 
plot was found for TTM3F, but with lower frequencies. 


D. LO-TO splitting vs temperature 

The frequencies of the librational and stretching modes 
are shown in table III HI Once again we compare our re¬ 
sults to experimental data^ii^i^S The comparison is im¬ 
perfect since the TIP4P/2005f and TTM3F data comes 
from data at finite k (the smallest k in the system). 
For all three systems (TIP4P/2005f, TTM3F, and ex¬ 
periment) the increase in the LO-TO splitting of the li¬ 
brational band is puzzling, since the right hand side of 
the LST relation predicts a decrease in splitting, cor¬ 
responding to a smaller dielectric constant and weaker 
dipole-dipole interactions. We found verifying the gener¬ 
alized LST equation is difficult because water contains ei¬ 
ther two or three Debye relaxations which must be taken 
into account.—!^ Uncertainties in how to fit the region 
of 1 - 300 cm“^ (.2- 9 THz), which includes contribu¬ 
tions from many H-bonding modes, precludes a direct 
application of the generalized LST relation to water. By 
ignoring this region, however, we were able to achieve an 
approximate validation of the generalized LST equation 
for TIP4P/2005f. A more detailed analysis of how to fit 
the low frequency region will be the focus of future work. 
Since the generalized LST equation is an exact sum rule 
it can be used to assist in testing the validity of different 
fit functions. 


E. Relation to phonons in ice 

Naturally we would like to find corresponding optical 
phonon modes in ice. As shown in figure [T] the dielec¬ 
tric spectra and LO-TO splitting of supercooled water 
resembles that of ice. Recently evidence has been pre¬ 
sented for propagating librational phonon modes in ice 
XL 42,43 Qf ^jj^g twelve librational modes of ice XI 
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FIG. 6. Imaginary parts of dielectric susceptibility. 

We compare (a) the non-polarizable model TIP4P/2005f, (b) 
the polarizable model TTM3F, and (c) experimental datai^ 
at 298 K. The effects of polarization can be seen in the LO- 
TO splitting of the stretching mode and in the low frequency 
featnres. 
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FIG. 7. Imaginary part of the distance decomposed 
susceptibility for TIP4P/2005f. Transverse (a) and lon¬ 
gitudinal (b) susceptibilities, calculated with a 4 nm box at 
300K, using the smallest k vector in the system. Gaussian 
smoothing was applied. Long range contributions to the li- 
brational peak extending to R — 2 nm are observed. 


are IR active (labeled WRl, RWl and RW2) and all three 
exhibit LO-TO splitting. The splittings have been found 
from Raman scattering to be 255, 135 and 35 cm“^42 
These modes all consist of coupled wagging and rocking 
motions. The WRl mode, which has the largest infrared 
intensity, most closely matches our results. WRl and 
RW2 have the same transverse frequency and RWl has 
a smaller infrared intensity, which may help explain why 
the librational band is well fit by a single optical mode. 
LO-TO splitting in the OH-stretching modes of ice Ih has 
been discussed previouslyJ^ 


F. Range of propagation 


The range of propagation of these modes can be calcu¬ 
lated as R = TVg where r is the lifetime and Vg = doj/dk 
is the group velocity. For TIP4P/2005f we find a range of 
propagation of ~ 1.1 nm for the longitudinal librational 
mode and ~ .3 nm for the transverse mode. Similar re¬ 
sults hold for TTM3F. 

To verify that the modes we observe are actually prop¬ 
agating and to further quantify the range of propagation 
we study the spatial extent of polarization dipole correla¬ 
tions as a function of frequency. We investigated several 
different methodologies that can be used to decompose 
a spectra into distance-dependent components (Supple¬ 
mentary Note 1). We choose to start with the polariza¬ 
tion correlation function: 

Hk,t) = (^Piik,0) ( 8 ) 

Here p^{k,t) represents either the longitudinal or trans¬ 
verse molecular polarization vector of molecule i. We 
now limit the molecules in the second sum to those in a 
sphere of radius R around each molecule v. 


<P{k,t,R) = ('^p,{k,0) ■ '^Pg{k,t) 


(9) 


j&Ri 


The resulting function exhibits the expected i? —?> 0 limit, 
yielding only the self contribution. R can be increased to 
the largest R in the system (•\/3L/2), where the full re¬ 
sponse function for the simulation box is recovered. As R 
increases, the contributions of the distinct term add con¬ 
structively and destructively to the self term, illustrating 
the contributions from molecules at different distances. 

Figure [ 7 ] shows the distance decomposed longitudinal 
and transverse susceptibilities for TIP4P/2005f in a 4 nm 
box at the smallest k available in the system. The entire 
region between 0 - 1000 cm“^ contains significant cancel¬ 
lation between the self and distinct parts, in qualitative 
agreement with a previous study4^ In the longitudinal 
susceptibility, the self component has two peaks (at 500 
and 900 cm“^ for TIP4P/2005f) representing the two IR 
active librational motions (rocking and wagging, respec¬ 
tively) . The self part is the same in both the longitudinal 
and transverse cases, reflecting an underlying isotropy 
which is only broken when dipole-dipole correlations are 
introduced. Further insight into the self-distinct cancella¬ 
tion comes from the results of Bopp, et. ah, who project 
the hydrogen currents into a local molecular frame, al¬ 
lowing them to study the cross correlations between the 
rocking and wagging librationsj^ They find that in the 
longitudinal case cross correlations between rocking and 
wagging contribute negatively in the region of 480 cm“^ 
and positively in the region of 740 cm“^, suppressing the 
lower frequency peak to zero and enhancing the higher 
frequency peak. 
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lesser extent. LO-TO splitting of about 100 cm“^ is ob¬ 
served in the 700 cm“^ librational peak. The results 
for acetonitrile are more ambiguous - we observe disper¬ 
sion in the broad peak at ~ 100 cm“^, however this peak 
contains contributions from translational and (free) rota¬ 
tional modes, as well as the CH 3 torsion mode, and it is 
not clear which modes are responsible for the dispersion 
(Supplementary Fig. 5). 


FIG. 8. Imaginary part of the distance decomposed 
longitudinal susceptibility for TTM3F at 300 K. Long 
range contributions are observed in the OH stretching band. 


In both the transverse and longitudinal cases as R in¬ 
creases a new peak emerges, corresponding to the prop¬ 
agating mode. Incidentally, the shift in the peak be¬ 
tween the self and distinct parts rules out the possibil¬ 
ity that the propagating mode is the proposed dipolar 
plasmon resonance since the dipolar plasmon must be 
a resonance of both the of single molecule and collec¬ 
tive motion4^“— Interestingly, there are very long range 
contributions to this peak. In our simulations with a 
4 nm box of TIP4P/2005f contributions persist up to 3 
nm in the longitudinal case and 2 nm in the transverse 
case. As noted, recent studies of ice XI suggest that the 
propagating modes consist of coupled wagging and rock¬ 
ing librations42ii^ The results for the transverse mode 
seem to confirm this hypothesis for liquid water, since the 
propagating mode peak lies between the single molecule 
rocking and wagging peaks. In the longitudinal case the 
propagating mode overlaps more with the wagging peak, 
suggesting a greater role for these type of librations in 
the longitudinal phonon. 


G. Methanol & acentonitrile 

To provide further evidence the aforementioned optical 
modes propagate through the hydrogen bond network of 
water we decided to repeat our analysis for other polar 
liquids, both H-bonding and non H-bonding. As an H- 
bonding liquid we choose methanol, which is known to 
contain winding H-bonded chains. According to results 
from MD simulation, most of these chains have around 
5-6 molecules^Si^, with a small percentage of chains con¬ 
taining 10-20 molecules4S Chain lifetimes have been esti¬ 
mated to be about .5 ps4S Therefore we expect methanol 
can also support a librational phonon mode that propa¬ 
gates along hydrogen bonds, but perhaps with a shorter 
lifetime and range than water. As a non H-bonding po¬ 
lar liquid we choose acetonitrile, because it has a struc¬ 
ture similar to methanol, but with the hydroxyl group 
replaced by a carbon atom. We find that the OH li¬ 
brational band of methanol (~ 700 cm“^— ) is indeed 
dispersive (Supplementary Fig. 4). As with water, the 
transverse spectra also exhibits dispersion, but to a much 


III. DISCUSSION 

In this work we have presented several lines of evidence 
for short lived optical phonons that propagate along the 
H-bond network of water. The longitudinal and trans¬ 
verse nonlocal susceptibility exhibit dispersive peaks with 
dispersion relations resembling optical phonons. As the 
temperature is lowered, the resonance frequencies and 
LO-TO splittings of these modes converge towards the 
values for phonons in ice Ih. By comparing our results 
with a recent study of ice XI we believe both modes likely 
consist of coupled wagging and rocking librations 

This work fundamentally changes our understanding 
of the librational band in the Raman spectra of water 
by assigning the lower and higher frequency peaks to 
transverse and longitudinal optical modes. Our analysis 
of the self-distinct cancellation indicates that the mid¬ 
dle Raman peak {k, 600 cm“^) belongs to the remnant 
of the single molecule wagging response which remains 
after the cancellation. We are also led to a new interpre¬ 
tation the librational region of the real part of the dielec¬ 
tric function. In the case of a lossless optical phonon the 
transverse phonon occurs where e'{uj) = 00 while the lon¬ 
gitudinal phonon occurs where e' {u) = 0. The presence 
of dampening smooths the divergence leading to a peak 
followed by a sharp dip. This is what is observed in the 
real part of the dielectric function of water between 300 
- 500 cm“^ (the features are shifted to lower frequencies 
by the tail of the low frequency Debye relaxation). 

One might wonder how our work relates to existing 
work on acoustic modes in water, in particular, the con¬ 
troversial “fast sound” modei^i^ Acoustic modes, which 
are observable through the dynamic structure factor, 
have been explored as means of understanding the hy¬ 
drogen bond structure and low temperature anomalies of 
water4 In this work we have argued that optical modes 
can also provide insight into water’s structure and dy¬ 
namics. The fast sound mode lies at much lower frequen¬ 
cies than the librational and OH stretch modes that we 
studied. The H-bond bending and stretching modes also 
primarily lie at at frequencies below the librational re¬ 
gion. However, normal mode analysis of liquid water and 
clusters shows that the H-bond stretching modes have a 
wide distribution of frequencies which overlaps with the 
librational modes, so some coupling between these modes 
is possibleRecently it was shown that there is cou¬ 
pling between the acoustic and optic modes in water - ie. 
between fluctuations in mass density and fluctuations in 
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charge densityi^ 

The large spatial range and coherent propagation of 
these modes is surprising and implies the existence of 
an extended hydrogen bond network, in contrast to ear¬ 
lier ideas about the structure of water which empha¬ 
size dynamics as being confined within small clustersi^ 
Simulations with larger simulation boxes are needed to 
fully quantify the extent of the longitudinal modes. The 
ability of water to transmit phonon modes may be rel¬ 
evant to biophysics, where such modes could lead to 
dynamical coupling between biomolecules, a phenomena 
which is currently only being considered at much lower 
frequenciesThe methodology used in this paper to 
analyse LO-TO splitting opens up a new avenue to un¬ 
derstanding the structure and dynamics of water. The 
fact that the librational LO-TO splitting increases with 
temperature instead of the expected decrease is likely 
due to significant changes in the structure of the liquid. 
One likely possibility is that the volume per “unit cell” 
term in equation [3] decreases with temperature. This 
could be caused by the local quasi-structure determined 
by H-bonding changing from a more ice-like structure 
(4 molecules per unit cell) to a more cubic structure (1 
molecule per unit cell). More research is needed to un¬ 
derstand the microscopic origin of the LO-TO splitting 
in water, both in the librational and stretching modes. 


IV. METHODS 

A. Theory of the nonlocal susceptibility 

If the external field is sufficiently small, the relation 
between the polarization response of a medium and the 
electric displacement field D for a spatially homogeneous 
system is given by : 

P{r,t)=eo f f dr'dt'X {'> — r',t—t')D{r',t') (10) 
Jv J-oo 

This expression Fourier transforms to: 

P{k,uj)=eox{k,uj)D{k,uj) (11) 

For isotropic systems, the tensor x can be decomposed 
into longitudinal and transverse components: 

X (k,uj) = XLik,uj)kk + XT{k,uj){I - kk) (12) 

The easiest starting point for deriving microscopic ex¬ 
pressions for XL{k,u}) and XT{k,uj) is the classical Kubo 
formula!^ 

R r°° d 

XL/T{k,u;) = ^ / dt-{PL/Tik,t) ■ P2/^(fc,0))e*“^‘ 
Co JQ ut 

(13) 

This expression relates the susceptibility to the time 
correlation function of the polarization in equilibrium. 


The longitudinal part of the polarization can be cal¬ 
culated by Fourier transforming the defining expression 
for the polarization V • P{r,t) = —p{r,t), leading to 
k ■ P = ip{k,t)/k = Pl- To calculate the transverse 
part of the polarization we use the method of Raineri & 
Friedman to find the polarization vector of each molecule 
(Supplementary Note 2)i^ We can rewrite eqn. [13] in 
terms of the normalized polarization correlation function 
(eqn. [7]), and taking into account the isotropy of water: 

poo 

XL/T{k,uj)=XL/T{.k,0) ^L/T{k,t)P‘^*dt (14) 

Jo 


B. Computational methods 

The three water models we used were TIP4P/e)S 
TIP4P/2005f,.^ and TTM3F.^ To simulate methanol 
and acetonitrile we used the General AMBER Force- 
field (GAFF))^ a forcefield with full intramolecular flex¬ 
ibility which has been shown to satisfactory reproduce 
key properties of both liquids.— Our TTM3F simula¬ 
tions were performed with an in-house code that uses 
the TTM3F force calculation routine of Fanourgakis and 
Xantheas. All other simulations were ran using the GRO- 
MAGS package (ver. 4.6.5).— We used particle-mesh 
Ewald summation for the long range electrostatics with 
a Coloumb cutoff of 2 nm for our 4+ nm simulations 
and a cutoff of 1.2 nm for our simulations with 512 
molecules. Our TTM3F simulations had 256 molecules 
and used Ewald summation with a Coulomb cutoff of .9 
nm. The principle TIP4P/2005f simulations contained 
512 molecules and were 8 ns long (Atout = 8 fs) and .6 
- 1.2 ns long (Afout = 4 fs). Other simulations were 1- 
2 ns long. Simulations with MeOH and ACN contained 
1,000 molecules and were 1 ns long. All simulations were 
equilibrated for at least 50 ps prior to outputting data. 

Because of periodic boundary conditions, the possi¬ 
ble k vectors are limited to the form k = 2TTnxi/Lx + 
2'KnyjILy 2Tmzk /where Ua,, nj,, and are inte¬ 
gers. We calculated correlation functions separately for 
each k and then average over the results for k vectors 
with the same magnitude, a process we found reduced 
random noise. 

One can question whether a purely classical treatment 
is justified here because the librational dynamics we are 
interested have frequencies of 700-900 cm“^ for which 
fiw ~ 3 — AksT at 300 K. Previously it was shown that 
the widely-used harmonic correction does not change the 
spectrum,Eurthermore, comparison of k resolved IR 
spectra taken from molecular dynamics and ab-inito DET 
simulation show that they give qualitatively similar re¬ 
sults for all frequencies below 800 cm“^— Eor the OH 
stretching peak, however, quantum effects are known to 
be very important. 


C. Fitting the librational band 

To obtain resonance frequencies and lifetimes for the 
librational peak in the imaginary part of the response we 
used a damped oscillator model. A Debye peak overlaps 
significantly with the librational band in both the lon¬ 
gitudinal and transverse cases and must be included in 
the peak fitting. Equation [14] can be used to relate the 
form of the time correlation function to the absorption 
peak lineshape. For Debye response one has the following 
expressions: 

Im{x(fc,w)} _ Alotd (15) 

X(fc,0) l + Uj'^T^ 

For resonant response with resonance frequency ujo{k) 
and dampening factor 7 = 1 /r we have: 

cos{ujot) 

lm{x{k^uj)} B / CJT LOT ^ 

x{k, 0) 2 \ 1 + (a; + 1 + (a; — a;o)^T2 y 

We find this lineshape (the Van Vleck-Weisskopf 
lineshape^i^) yields results identical to the standard 
damped harmonic oscillator response for the range of 
r, a;o values we are interested in. We found a two func¬ 
tion (Debye -h resonant) fit worked very well for fitting 
the librational peak in the longitudinal case (Supplemen¬ 
tary Figs. 7 and 8 ). The H-bond stretching peak at ^^200 
cm“^ overlaps with the librational band for 2 < k < 2.5, 
and we found that it can be included in the fit using an 
additional damped harmonic oscillator, but usually this 
was not necessary. Because of this overlap and due to the 
broad nature of the transverse band, the fitting in the 
transverse case is only approximate. We found this was 
especially true for TTM3F and the experimental data, so 
we do not report lifetimes for such cases. 
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VI. SUPPLEMENTARY INFORMATION 



t(ps) t(ps) 


FIG. 1. Longitudinal polarization relaxation functions. Shown for 512 TIP4P2005/f (left) and TTM3F (right) at 300 

K. 
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Supplementary Table 1: Experimental peaks in Raman, dielectric, and IR spectra. This table shows the 
correspondence between 3 peak Raman fits and 2 peak dielectric/IR fits to the librational region at 298 K. 





















FIG. 2. Longitudinal (left) and transverse (right) relaxation times for 512 TIP4P/2005f. Computed for the 
underlying exponential of the relaxation. The points are interpolated by Akima splines. The transverse relaxation time here is 
equal to the Debye relaxation time, which at A: = 0 is 11 ps at 300 K for TIP4P/2005f. Experimentally it is 8.5 ps.— 


A. Supplementary info: Dispersion relations and dampening factors 
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FIG. 3. Longitudinal (left) and transverse (right) dispersion relations (circles) and dampening factors (squares) 
for 512 TIP4P/2005f. These curves were obtained from a two peak (Debye + resonant) fit. In contrast to the longitudinal 
mode, the transverse mode is much more damped. However, the dampening factor changes significantly with temperature, 
also in contrast to the longitudinal case, and at 250 K becomes relatively small. Beyond 2 A the peak due to the damped 
resonance starts to disappear so values beyond 3 A are not shown. 
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VII. SUPPLEMENTARY INFO: METHANOL &. ACETONITRILE 



FIG. 4. Longitudinal (left) and transverse (right) dielectric susceptibility for a system of 1,000 MeOH molecules. 

The longitudinal librational peak at « 700 cm“^ clearly disperses with k, while the transverse peak at « 600 disperses 

slightly with k. The higher frequency peaks exhibit no dispersion. The static dielectric function e{k, 0) has not converged 
properly in the transverse case, so the magnitude of the peaks is not converged. 


CO (THz) “ (THz) 



FIG. 5. Longitudinal (left) and transverse (right) dielectric susceptibility for a system of 1,000 acetonitrile 
molecules. The broad band which peaks at 100 cm“^ exhibits dispersion. We hypothesize this dispersion is due entirely to 
the translational modes, however we cannot say for sure since the librational and translational modes overlap in this region. 
The peak at « 500 cm“^ is due to GCN bending. The static dielectric function e(fe, 0) has not converged properly, so the 
magnitude of the transverse peaks is not converged correctly, but the position of the peaks and dispersion can be seen. 
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A. Supplementary info: Examples of fitting 



FIG. 6. Example fits of the transverse susceptibility of TIP4P/2005f at 300 K. Fit with a Debye function and one 
damped harmonic oscillator at fc = .25A ^ and k = l.lA The residual show the parts not captured by the fit. 




FIG. 7. Example fits of the longitudinal susceptibility of TIP4P/2005f at 300K. Fit with a Debye function and one 
damped harmonic oscillator at fc = .25A ^ (left) and k = l.lA ^ (right). Two peaks appear in the residual - the lower frequency 
peak is dispersive, having the same dispersion relation as the fitted peak, suggesting that it is actually part of the dispersive 
peak lineshape that is not captured by our lineshape function. The higher frequency peak in the residual is non-dispersive and 
is in the same location for both the transverse and longitudinal susceptibility. 
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B. Supplementary info: spatial decomposition of spectra 


There are several different ways to decompose a spectra into contributions from molecules separated by distance R\ 

Kirkwood dipole-sphere method 

This is the method we choose, which is a modification of the sphere-sphere method (see below). We start with the 
time-correlation function of interest : 

( 16 ) 


Here fi can be replaced with any dynamical variable of interest, for instance t) or j{t). We omit the k dependence 
for simplicity. 

The most straightforward way is to limit the molecules around each molecule i to those in a sphere of radius R: 


(j){t,R) 



(17) 


This is similar to the method employed by Bopp & Kornyshev. During the the course of a simulation molecules 
enter and leave each sphere, which creates noise, requiring longer averaging times. This can be improved by utilizing 
a smooth cutoff function: 


where 






1 

“ 1 + g(R.i-R)/D 


(18) 


(19) 


Here D is a sharpness parameter determining the relative sharpness of the cutoff. We choose not to use smoothing 
however, finding it to be unnecessary. The result is a spectra x(A;,w,i?) showing contributions from molecules up to 
radius R. The resulting function exhibits the expected i? —>■ 0 limit, yielding only the self contribution. In the i? —>■ oo 
limit, the original full response function is recovered. This function can then be numerically differentiated to show 
the contributions from shells of thickness Ai? centered at distance R. 


Sphere-sphere method 

Another method discussed by Heyden, et. al. (2010) is to study the autocorrelation of the total dipole moment of a 
sphere of radius R centered around a reference molecule, and then average over each molecule in the system.— 

= ( 6 ) ■ (^)) ( 20 ) 
i 

where 

(t) = M(t) X (21) 

j€Ri 

Heyden, et. al. recommend the normalization factor Mi{t) = (1 -I- lo normalize for number of 

molecules in each sphere. This normalization factor is chosen so that in the bulk limit (i? —>■ oo) the original full 
response function is obtained (in that limit A/) = ^/y/Nrao\)- In limit i? —^ 0 only the self-term contributes. 
Results from this method must be interpreted with a bit of care since the calculation includes all cross-correlations 
between molecules within the sphere centered around the reference molecule. We found that this method is more 
sensitive to intermolecular correlations, in particular the H-bond stretching at ~ 250 cm“^ (not shown). Altogether 
though we found the results from this method are complementary with our results from the dipole-sphere method. 

Spatial grid method 

To achieve higher resolution, Heyden, et. al. also introduce a spatial grid method;^ The method works by bining 
the molecular dipoles into grid cells. To reduce noise caused by moleules moving in and out of bins the binning is 
Gaussian, meaning the dipoles are smeared with a Gaussian function. Unlike the other methods the spatial grid 
method does not yield the self part as i? —>■ 0 so this limit requires special interpretation. 
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C. Supplementary info: calculation of polarization vectors 


Bopp & Kornyshev show that to get accurate results in k space it is important to use the polarization vectors for 
each molecule rather than just the dipole moment. To calculate the polarization vector we use the method of Raineri 
& Friedmani^ We utilize the defining relation for the polarization: 

V-P(r,t) = -p(r,t) (22) 

When transformed into Fourier space this becomes: 

ik ■ P{r,t) =—p{k,t) (23) 

We introduce polarization vectors for each molecule Pi{k) so that we have 

(24) 

i 

where 

zfc.p,(fe) = -^g„e-*"- (25) 


The molecules are indexed by i and the atoms on each molecule are indexed by a. r^a = ‘fi(t) — fait) is the distance 
from each atomic site to the center of mass of molecule i. Following Raineri & Friedman, we use the identity 


= 1 + a: 


dse^ 


and taking into account the charge neutrality of each molecule we obtain 


Pi(fc) = - 


dse 


-ih-VaiS 


(26) 


(27) 


Pi(fe) = 


QaTc 


^ ik-Vr 




(28) 


The transverse part is then calculated as Pt = k x P, while the longitudinal component is Pl 
longitudinal component can also be calculated more directly from: 




This yields the following Kubo formula for the longitudinal part of the response: 


B d 

XL{k,uj) = dt—{p{k,t)p*{k,0))e^‘^* 

For a system composed of point charges, the charge density is : 

P{r, ^) = -^ XIX! 

i OL 


k ■ P. The 


(29) 


(30) 


(31) 


Again, the index i runs over the molecules while a runs over the atomic sites on each molecule. The charge density 
in k-space becomes: 


p{k,t) 


V 






Note that this can be Taylor expanded as: 


p{k,t) 

k 


1 {—ik ■ ria{t))‘^ 

k nl 

i a n 

M{k,t) + Q{k,t) + 0{k,t) + ■ ■ ■ 


(32) 


(33) 
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FIG. 8. Imaginary part of the longitudinal (top) and transverse (bottom) polar structure factor for TIP4P/2005f at 250 K, 300 
K, 350 K, and 400 K (left to right). Note the increased intensity of the low frequency, high wavenumber intramolecular mode 
at higher temperatures. This is likely due to weaker H-bonding and greater freedom for inertial motion, which is responsible 
for this band. 


Here M{k, t), Q{k, t), 0{k, t) are contributions due to the molecular dipoles, quadrupoles and octupoles. In the limit 
/c —^ 0 it from supplementary equation [30] it can be seen that only the dipole term contributes to the susceptibility. 
In the fc —>■ 0 limit one obtains 

dt-{ML{k,t)-Ml{k,0))e^^^ (34) 

with 

'^mol 

Mi(fc,t)= ^ (35) 

This type of expression has been used previously as an approximate expression at small However, Bopp & 
Kornyshev show quite convincingly that for water the higher order multipole terms are very important, even at the 
smallest k available in computer simulation.— Neglect of the higher order terms leads to severe consequences at large 
fc, and one will not recover the physical limit lim XL/rik) = 1 unless higher order terms are included. 

fe—>oo ' 

D. Supplementary info: polarization-polarization structure factors 

To assist in visualizing XL,Tik,u}) we introduce longitudinal and transverse “polarization-polarization structure 
factors: 

nOO 

S[^ik,uj)= / ^L,T{k,t)e^"'dt (36) 

do 

Thus, XL,Tik,uj) = XL,T{k,0)S^^{k,uj). These plots are shown solely because they provide a nice visual overview of 
the features in the nonlocal susceptibility. The main novel feature that appears in these plots is the low frequency 
acoustic-like mode originating at ~ 60 cm“^. This mode is purely intramolecular in nature and arises from inertial 
rotation i^^^^ At very high wavenumbers {k > 7)) the relaxation is described by a rapidly decaying exponential and a 
Gaussian function— 

$L(fc, t) = A(fc)e*/^i('=) -k (37) 

Gaussian relaxation yields the following equation for the imaginary part of the susceptibility: 

$(fc,t) = 

Im{x(fc, w)} = x(fc, 

The real part is: 

Re{x(fc, w)} = xik, 0)H(t - t‘^ujF{ujt/2)) 


(38) 

(39) 
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where F{) is Dawson’s integral. The Gaussian form for the correlation function can be derived by considering a free 
rigid dipole subjected to Brownian kicks. In that case the relaxation function can be computed exactlyi^ 


= exp 


t 

n 


II 

n 


11 — exp 



(40) 


where ti = ^/2kBT and T 2 = I and ^ is the friction. In either the limit ^ 0 or t —>■ 0 one obtains the Gaussian 

form. Thus the interpretation of the Gaussian form is that it is due to fast inertial relaxation. It has been suggested 
that such inertial relaxation is origin of the Poley absorption that has been found in some dipolar liquids.— “Poley 
absorption” appears to be used as a general term for absorption of unknown origin found in many polar liquids around 
10 cm“^ (.3 THz))^ first described by Poley in 1955— It has been variously described as being due to fast inertial 
“rattling” of molecules within their potential energy basins or as due to fast librational/inertial motion analogous to 
the rotational absorption of gas molecules— 




